Nonadditive hard-sphere fluid mixtures: A simple analytical theory 



O 
(N 

o 

o 



o 



I 

O 

o 



> 

m 

o 

^, 

l> 

o 



X 



Riccardo FantonO 
National Institute for Theoretical Physics (NITheP) and Institute of Theoretical Physics, Stellenbosch 7600, South Africa 

Andres Santof[J 
Departamento de Fisica, Universidad de Extremadura, E-06071 Badajoz, Spain 

(Dated: October 14, 2011) 

We construct a non-perturbative fully analytical approximation for the thermodynamics and 
the structure of non-additive hard-sphere fluid mixtures. The method essentially lies in a heuristic 
extension of the Percus-Yevick solution for additive hard spheres. Extensive comparison with Monte 
Carlo simulation data shows a generally good agreement, especially in the case of like- like radial 
distribution functions. 
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I. INTRODUCTION 

The van der Waals ideas [l| show that the most im- 
portant feature of the pair potential between atoms or 
molecules is the harsh repulsion that appears at short 
range and has its origin in the overlap of the outer elec- 
tron shells. These ideas form the basis of the very suc- 
cessful perturbation theories of the liquid state. This, 
along with fruitful applications to soft matter [2|, ex- 
plains the continued interest in hard-sphere reference sys- 
tems 0. 

The simplest model for a fluid mixture is a system of 
additive hard spheres (AHSs) for which the like-unlike 
collision diameter (aij) between a particle of species i 
and one of species j is equal to the arithmetic mean 
af^'^ = \{(Tii -\- CTjj). A more general model consists 
of nonadditive hard spheres (NAHSs), where the like- 
unlike collision diameter differs from crff'^ by a quantity 
Aij = {aij — crf^'^) / (jf^'^ called the nonadditivity param- 
eter. As mentioned in the paper by Ballone et al. [j], 
where the relevant references may be found, experimental 
work on alloys, aqueous electrolyte solutions, and molten 
salts suggests that homocoordination and heterocoordi- 
nation [3, 0] may be interpreted in terms of excluded vol- 
ume effects due to nonadditivity (positive and negative, 
respectively) of the repulsive part of the intermolecular 
potential. NAHS systems are also useful models to de- 
scribe real physical systems as rare gas mixtures [7[ and 
colloids [a-fllj. For a short review of the literature on 
NAHSs up to 2005 the reader is referred to Ref. ^. 

The well-known Percus-Yevick (FY) integral-equation 
theory [l[ is exactly solvable for a mixture of three- 
dimensional (3D) AHS mixtures [Tj, [ij]. The solution 
has been recently extended to any odd dimensionality 
[iMi- On the other hand, any amount of nonadditivity 
(Aij 7^ 0) suffices to destroy the analytical character of 
the solution and so one needs to resort to numerical meth- 



ods to solve the FY or other integral equations Q- 

The aim of the present paper is to propose a non- 
perturbative and fully analytical approach for 3D NAHS 
fluid mixtures, which can be seen as a naive heuristic 
extension of the FY solution for AHS mixtures. In do- 
ing this, we are guided by the exact solution of the one- 
dimensional (ID) NAHS model [16l - [l9J and some physical 
constraints are imposed: the radial distribution function 
(RDF) gij(r) must be zero within the diameter aij, the 
isothermal compressibility must be finite, and the zero 
density limit of the RDF must be satisfied. We find that 
this strategy gives very good results both for the thermo- 
dynamics and the structure, provided that some geomet- 
rical constraints on the diameters and the nonadditivity 
parameter are satisfied. This makes our approach partic- 
ularly appealing as a reference approximation for integral 
equation theories and perturbation theories of fluids. 

The paper is organized as follows: in Sec. |ll] we de- 
scribe the NAHS model outlining the physical constraints 
that we want to embody in our approach. The latter is 
constructed by a three-stage procedure (approximations 
RFA, RFA+, and RFA^^"^) in Sec. HhI In Sec. |IV] we 
present the results for the equation of state from our ap- 
proximation, comparing them with available Monte Carlo 
(MC) simulations. The results for the structural proper- 
ties are presented in Sec. |Vl where we compare with our 
own MC simulations. Finally, Sec. lVIl is devoted to some 
concluding remarks. 



II. THE NAHS MODEL 

An n-component mixture of NAHSs in the d- 
dimensional Euclidean space is a fluid of A'^; particles of 
species i (with i = 1,2, ... ,n), such that there are a to- 
tal number of particles N = X]i=i ^i in a volume V, and 
the pair potential between a particle of species i and a 
particle of species j separated by a distance r is given by 
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pairs (i, j) we recover the AHS system. In a binary mix- 
ture {n = 2), Ai2 = A21 = A is the only nonadditiv- 
ity parameter. If A = —1 one recovers the case of two 
independent one-component hard-sphere (HS) systems. 
In the other extreme case cti = CT2 = with (T12 finite 
(so that A —J- 00) one obtains the well known Widom- 
Rowlinson (WR) model [20;l21|. Another interesting case 
is the Asakura-Oosawa model [22], [2j] (where 172 = and 
A > 0), often used to discuss polymer colloid mixtures 
and where the notion of a depletion potential was intro- 
duced. The NAHS system undergoes a demixing phase 
transition for positive nonadditivity [24l - l28l | . A demixing 
transition might also be possible, even for negative non- 
additivity [23. ISOj. provided the asymmetry ratio oxjai 
is sufficiently far from unity. In the present paper we will 
only consider the NAHS system in its single fluid phase. 
Let the number density of the mixture be p = N jV 
and the mole fraction of species i be Xi = Pi/p, where 
Pi = Ni/V is the number density of species i. From these 
quantities one can define the (nominal) packing fraction 
rj = VdpMd, where Vd = {■K/4)'^/^T{l + d/2) is the volume 
of a d-dimensional sphere of unit diameter and 



Mk 



E-- 



(2.2) 



denotes the fcth moment of the diameter distribution. 
The NAHS model, in the thermodynamic limit N — >■ 

00 with p = N/V constant, admits an analytical exact 
solution for the structure and the thermodynamics in d = 

1 [16l - [l9l |. Moreover, the AHS model in odd dimensions 
is analytically solvable in the PY approximation [l3l - [l5l | , 
the result reducing to the exact solution of the problem 
for d = 1 but not for d> 3. 



A. Basic physical constraints on the structure 

The RDF gij{r) must comply with three basic condi- 
tions: 

(i) gij{r) must vanish for r < Uij. More specifically, 
for distances near an, 



gtj{r) = e(r - (j,j) [g,j(cr,+ ) + g[j{(T^){r - a,j) + 
where Q(x) is the Heaviside step function. 



(2.3) 



(ii) In the fiuid phase the isothermal compressibility 
X must be finite. This implies (see below) that 
the Fourier transform hij (q) of the total correlation 
function hij{r) = gij{r) — 1 has to remain finite at 
q = or, equivalently, 

/•oo 

/ dr r^h.j (r) = finite for < a < d - 1. (2.4) 

Jo 



(iii) In the low density limit, the RDF is 

lim g,,(r) - g-^-''-'/'^^^ = e(r - a„), (2.5) 

ks and T being the Boltzmann constant and the 
absolute temperature, respectively. 

As a complement to Eq. (|2.5I) . we give below the exact 



expression of gij (r) to first order in density 31 



fc=i 
X (r - o-jfe - akj)"^ [r^ + 2{aik + <Jkj)r 
-3(f7.fc-afcj)2]+0(p2)|. (2.6) 

B. The two routes to thermodynamics 

For an athermal fiuid like NAHSs there are two main 
routes that lead from the knowledge of the structure to 
the equation of state (EOS) [l|. These may give different 
results for an approximate RDF. 

The virial route to the EOS of the NAHS mixture re- 
quires the knowledge of the contact values gij [af, ) of the 
RDF, 

od-l " 
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where Z = p/ pksT is the compressibility factor of the 
mixture, p being the pressure. 

The isothermal compressibility %, in a mixture, is in 
general given by 



1 [dp 






ksT \dpjj,^,^^y ksTf-^ \dpjj.^^^y 



1 - p X XiXjCij{0), 



(2.8) 



where Cij (q) is the Fourier transform of the direct corre- 
lation function Cij{r), which is defined by the Ornstcin- 
Zernike (OZ) equation 

n 

htj{q) == c^J{q) + X Pkhik{q)ckj{q)- (2.9) 

fe=i 



Introducing the quantities hij{q) = y/piPjhij{q) and 
Cij(q) = ^piPjCij{q), the OZ relation becomes, in ma- 
trix notation, 



c(q)-h(<7). l + h(9) , 



(2.10) 



where I is the n x n identity matrix. Thus Eq. (I2.8P can 
be rewritten as 



"=E/^ 



I'^j I'-'tj ^tj 
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l + h(0) 



(2.11) 



In Eq. (|2.11l) . and henceforth, we use the notation A^^ 
to denote the ij element of the inverse A~^ of a given 
square matrix A. 

In the particular case of binary mixtures {n — 2), Eq. 
(P1T|) yields 

_ [1 + pXihii{Q)][l + /9X2fe22(0)] - p^XiX2feJ2(0) 

1 + pa:iX2 [/ill (0) + 7^22 (0) - 2/^12(0)] 

(2.12) 
The compressibility route to the EOS can be obtained 
from 



Z^(77)= / dxx-\vx). 



(2.13) 



C. The one-dimensional system 

The exact solution for nonadditive hard rods [d = 1) 
is known ^19, 32„i3j|- First, let us introduce the Laplace 
transform 

/•OO 

G,,{s)^ / dre-^^g,,{r). (2.14) 

In terms of this quantity the exact solution has the form 

_ (2.15) 

where 



n 

G,,(s)-^=^P.fc(s)g,,(s) 

V^«^J k=i 



Pij{s) = ^XiXjK.j- — - — (2.16) 

is proportional to the Laplace transform of the nearest- 
neighbor probability distribution and 



Q{s)^[\-pP{s)\ 



(2.17) 



In Eqs. ([2T5)) and ((2J6)) . ^ ee p/keT = pZ, while 
Kij — Kji are state-dependent parameters that are de- 
termined as functions of ^ from the condition (j2.4p , which 
implies linis_+o sGij{s) ~ 1, as well as requiring the ratio 
Kij/Kik to be independent of i [1^. Those conditions 
also provide the exact EOS in implicit form, i.e., p as a 
function of ^. 

Of course, the above results also hold for additive hard 
rods. In that case, the additive property o-y = af'^'^ = 



2 {(Ji+aj) allows us to rewrite the solution in other equiv- 
alent ways. To that end, let us define 



so that 



U, = K,,e-i< 



ij V'^J — ■y^^i'^j-'-^i^ 



s + e ' 



where 



Q-.\.) = e-^yfJ^a,(.), 



0"i3 = 2*-^' ^'^•''^ 



a,(.)^(l + ^)5.,-^L.,e--^ 



Here we have made use of the property 



„ .-..add I ^ 



and 



It is easy to prove that 



),,(.) = e-^%/|^C^^(.), 



(2.18) 



(2.19) 



(2.20) 



(2.21) 



(2.22) 



(2.23) 



(2.24) 



thanks to the property aik + auj = Uij . Consequently, in 
the additive case, Eq. (|2.15p becomes 



Gtri-s) 



^L,feC,-/(s), (2.25) 



fe=i 



where use has been made of the additivity property 
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flfcj = 0-. 
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The additive solution turns out to be 

^ 1 



u 



p 1 - pMi 



(2.26) 



(2.27) 



The fact that Lij = const allows one to rewrite Eq. (|2.25p 
in yet another equivalent form, 



Gr(^) = ^^E^-'^-^fc.'(^)' 



(2.28) 
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where 



Bij{s) = Sij -L,jipQ{a,s) 



Cij{s) + - {xi -%) . 



In the first equality. 



ipQ{x) = e ''-I. 



(2.29) 
(2.30) 



While lims_i.o sCij{s) — ^ {Stj — Xi) ^ 0, but det[sC(s)] — 
0{s), in the ease of the matrix B(s) one has 
Imig^o sBij(s) — 0. On the other hand, in both 
cases, lmis^ooCij{s) = linis-foo Bij{s) ~ Sij, so that 

Hms^oose'^'j Gij{s) = Ljj = ^ /p. 

It turns out that Eqs. ([JT^ . 1^^ . and 1^^ are also 
obtained from the PY solution for additive hard rods. 
Thus, the PY equation yields the exact solution in the 
additive case, but not in the nonadditive one. 

It is important to bear in mind that, if one inverts 
the steps, it is possible to formally get Eq. (|2.15|) from 
Eq. p.25p . In other words, starting from the form ()2.25p 
of the PY solution for the ID AHS system, allowing Lij 
and ^ to be free, and carrying out some formal manipula- 
tions, one arrives at an equivalent form, Eq. (j2.15|) . that, 
if heuristically extended to the NAHS case {dij ^ crfj^*^), 
coincides with the exact solution to the problem. How- 
ever, it is not possible to recover (|2.15p starting from the 
form (|2.28p since the property Ly = const, only valid in 
the additive case, cannot be reversed. 



where L(s) and B(s) are matrices given by 

L,,{s)=L^Ul'^'>s, 



(2.37) 



Btj (s) = S,j 
where 



2'KpXi 



L[fip2{(Tis) + L^fsipiiais) 



tpi{x) ^ e ^-l + a:, (p2{x) 



(2.38) 



1 + x-—. (2.39) 



Similarly to the ID case, linis^o sBij(s) = 0. In fact, 
Eqs. ([236)) - ((238)) are the 3D analogs of Eqs. ([2:28)) and 
(|2.29|) . For the general structure of the PY solution with 
d ~ odd, the reader is referred to Ref . [1^ . 

Also as in the ID case, limg-j-oo Bij (s) = Sij and so, 
according to Eq. (P3^ . 
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(2.40) 



D. PY solution for three-dimensional AHSs 

In this subsection we recall the PY solution for AHSs 
in three dimensions {d ~ 3) p^. [l4|. 

First, one introduces the Laplace transform of rgij{r), 

G,j{s)= / dre-^Vg,,(r). (2.31) 



From Eq. ((23)) it follows that 

se^^^'Cjis) = fT„-,9,j(4) + [gy(4) + fT,J5^(4)] s^^ 
+0{s-^). (2.32) 

Next, Eq. ()2.4p implies, for small s. 



.^G,,(.) = l + 4').2 + i/«.3 + . 



(2.33) 



with H\f = finite and H\f = -/ijj(0)/47r = finite, 
where in general 



H. 



(a) 



dr{-rYrhij{r). (2.34) 



Further, in view of Eq. ()2.33p . the coefficients of s^ and s 
in the power series expansion of s'^Gij{s) must be 1 and 
0, respectively. This yields 2n^ conditions that allow us 
to find 111 



(0) 



01 + 6*2(7 
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720'j(T_,-, 



(2.41) 



where 9i = 1/(1 - t]) and 6*2 = ^{Kh j Mz)^ j {\ - t])^. It 
is straightforward to check that Eq. p.36p complies with 
the limit (|05)) . 

The expressions ()2.7p and (|2.13p which follow from the 
solution of the PY equation of AHS mixtures are 
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(2.43) 
Usually, the virial route underestimates the exact results, 
while the compressibility route overestimates them. 



Finally, Eq. ([23]) yields 



Um G.ij{s) 






(l + CTys). (2.35) 



Equations (|23T)) - (|235)) hold both for NAHSs and AHSs. 
The PY solution for AHSs can then be written as [IJ, 



El 



GTi^ 



■Y,L,k{s)B-^{s), (2.36) 



k=l 



III. CONSTRUCTION OF THE 
APPROXIMATIONS 

As stated in Sec. HI the main aim of this paper is to 
construct analytical approximations for the structure and 
thermodynamics of 3D NAHSs. On the one hand, the 
approximations will be inspired by the exact solution in 
the ID case (see Sec. Ill C)) . On the other hand, they will 
reduce to the AHS PY solution (see Sec. IIIDj) . More- 
over, as a guide in the construction of the approximations 
and also to determine the parameters, the basic physical 



requirements ((231) - (|23|) [or, equivalently, ((232)) . ((2?33l) . 
and (|2.35p ] will be enforced. 

The driving idea is to rewrite Eq. ()2.36|) in a form akin 
to that of Eq. (|2.15|) . by inverting the procedure followed 
in Sec. Ill CI This method faces several difficulties. One 
of them is that, as said before, Eq. ()2.36p is the 3D analog 
of Eq. dOSl), but not of Eq. (g^Sl), and it is not possible 
to recover directly (i.e., without further assumptions) Eq. 
(I2T5)) from Eq. ((2^ . One could first try to rewrite Eq. 
(|2.36p in a form akin to that of Eq. (|2.25p . i.e., a form 
where the matrix B given by Eq. (|2.38p is replaced by a 
matrix C such that Imis^Q sCij (s) 7^ 0. But, given the 
intricate structure of Eq. (|2.38p and the fact that neither 

L}/ nor L}/ are constant, this does not seem to be an 

''J ^j ^^^_^ 

easy task at all. Therefore, we will work from Eq. (|2.36p 
directly. 



The AHS PY solution revisited 



First, define 



Pij{s) EE y^xlxje "'^ ''Lij{s), 



Q^j{s) 



-Br^is), 



(3.1) 



(3.2) 



so that 



Q^\s) = e" 



'-B,,[s) 



L^!f •^2{(Tis) + L')^hipi{ais) 



(3.3) 



Inserting Eqs. ((3T|) and ((3?2]) into Eq. ((236)) we finally 
get 



G,,{s) 



V^^ k=i 



y2p.k{s)Qk,{s), 



(3.4) 



where use has been made of the additive property ()2.26p . 
We emphasize that Eq. p.4p is fully equivalent to Eq. 
p.36p and thus it represents an alternative way of writing 
the PY solution for AHSs. In both representations the 
coefficients L^ ■ and L,- are given by Eq. (|2.4ip . On the 



other hand, since the structure of Eq. p.4p is formally 
similar to that of the exact solution for ID NAHSs, Eq. 
(|2.15p . it might be expected that Eq. (|3.4p is a reasonable 
starting point for an extension to 3D NAHSs. 



B. Approximation RFA 

1. The proposal 

A possible proposal for the structural properties of 
NAHSs is defined by Eq. ((33) with 



ij V*^/ — -s/^i^j^ ^ij\^)-) 



(3.5) 



Q^■,\s) = 6,, + 



2ttp^/x-x- ^^ 



Lf^^2{h,s) + Lf^s^i{h,s) , (3.6) 



where Lij{s) is still given by Eq. (|2.37p [with i,- and 



L\,' yet to be determined] and 



Oij — (7ij + dij . 



(3.7) 



Equations (|3.5p and p.6p are obtained from Eqs. p.ip 
and (|3.3p . respectively, by the extensions af^'^ —>■ <Jij and 
(Ji -^ bij [compare Eqs. ((^^ and ((XT)) ] . Note that Eq. 
(13.61) can also be written as 



Qi^i^ = h - 2^P^^ [iV,,(,s)e''-- - i«,(.s)e-'^-' 



where 



(3.8) 



= r(o) 



N,,{s)^l4^>\l-h,s^ 



I^s (1 - h,s) . 



(3.9) 
Of course, the coefficients L^- and Li-' are no longer 
given by Eq. (|2.4ip but are obtained from the physical 
conditions 



lira s Gij{s) = 1, 



s->-0 



lims-i [s'^Gij{s)^l] =0, 



s-J-O 



(3.10) 



(3.11) 



which follow from Eq. (|2.33p . To that purpose, it is con- 
venient to rewrite Eq. p.4p as 



s^ ^ \/xiXkGtk{s)Qk^{s) = Pij{s). 

k=l 

Using Eqs. ((S3|) and ((SD), Eq. ((XTU)) implies 



k=l 



(3.12) 



(3.13) 



Likewise, Eq. p. lip gives 



'"^^kbl 



TTp2_^Xkblj 

fe=i 



akj L 



(1) 



lJ^}bk. 



^..[l^-\l^K 



= Lif-a.,Ll°'.(3.14) 



Equations ((3.13P and p.l4p imply that both L\,j and 
-(1) ^..r(0) 



L-- — aijLi- are independent of the subscript i, i.e. 



L 



ij — '-'3 ' ^ij — ^j + '^*j ^j ' 



(0) 



(3.15) 



where Sj and Tj are determined from Eqs. ()3.13p and 
((XTI)) . The solution is 



5", = 7- 
' 1 



1 - np'^j 



(1 - TTpAj) (1 - 7rp*j) - n^ p'^ pj\2,o^j 



(3.16) 



T,= 



Trpflj 



^ (1 - TrpAj) (1 - 7rp*j) - 7r2p2^j.|2 pf7j 

where we have called 

, _ 1 



*J = 3^*^13,0 -Mj|2,l' 



and 



^j — A'j|3,l 



f^j\p,i 



1 



f^j\2.2 



7^^j\i.,o, 



Y.'^h^lAr 



k=l 



In the additive case {hkj ~ Ufe) one has A^ = g -11:13 



(3.17) 



(3.18) 



(3.19) 



(3.20) 



(3.21) 



M, 



^Msa^, so 



\M2aj, *j = iMa - \M2Gj, and rj^ 

that Sj = 6'i +02(Tj and T^- = ~5^20'|, in agreement with 
Eq. (|2.4ip . In the case of binary nonadditive mixtures 
(A 7^ 0), it can be easily checked that the common de- 
nominator in Eqs. p.l6p and (|3.17l) is positive definite. It 



only vanishes if A = — 2(T2/(ci + 0-2) (assuming (T2 < ^i) 
and 77 = 1 + X2cr2/xiaf. 

Equation ()3.15p closes the approximation p.4p - p.6p . 
It relies on the same philosophy as the so-called rational- 
function approximation used in the past for HS and re- 
lated systems [Tsl, [SJ] and, therefore, we will use the 
acronym RFA to refer to it. The explicit forms of Gij{s) 
for binary mixtures (n = 2) are presented in Appendix 

El 



2. Low-density behavior 
To first order in density, Eqs. ((3T5l) - ((3Tf|) yield 



Lf>^^l + npA,+0{p'), 



(3.22) 



4f = a„- + np (a.jAj + Q/) + O(p^). (3.23) 



Thus, 



2TTp^/XiXj 

)y(s) = Sij g e"'^" [f2{bijs) + aijS(pi{bijs)] 



+0{p'). 
Insertion into Eq. p.4p yields 



(3.24) 



G.J is) 



[1 + a,,s) + ^p'-^ [A, + (a., A, + ^,) -s] - ^ ^ Xfee-('^-+'^'=^> 

fc=i 



2Tip 



X (1 + a,fcs)(l + (Tk,s) + ^Y1 2:fce-('^--'"=^)^(l -t- a.fes) 
+0{p'). 



k=l 



l-akjs- -{alj -alj] 



Laplace inversion gives 

n 

9tjir) = e(r - cr.y) + ^e(r - cr,j) (Aj-r + ^j) - T:r-^^XkQ{r - a,k - crfcj)(r - cr,,fe - crfe^)^ 



fc=l 



X P -I- 2(crife -I- (TfcjOr - 3{aik - o-fcj)^] -f — ^ a:fce(r - CTife -I- akj)ir - aik + Ukj) 

k=l 

X [r^ + {(7.ik - akj)r^ - {^crji. + Gal^ + 2aikakj - alj)r 
+3(a,, + ak.Kaf, + a% - 2a%)\ + Oip"). 



r 



(3.25) 



(3.26) 



As a consequence, approximation RFA is consistent 
with the exact limits p.Sp and (|2.35p . To first order in 
density, the approximation correctly accounts for singu- 
larities of gij (r) at distances r = <Tij and r = aik + (^kj , 
k = 1, . . . ,n [see Eq. p.6p ]. On the other hand, we see 
from Eq. p.25p that, already to first order in density, 
approximation RFA introduces spurious singularities at 
r = Oik ^ cikj 7^ cTy- Oug might even have dij^k = 



Oik—akj—Oij < 0. In particular, da-^k — ofk Aik becomes 
negative if Aik < 0. Analogously, dij. 



^add 



•y;'j 



Aij IS 



negative if A^ > 0. Therefore, approximation RFA does 
not verify in general the condition (|2.3p . It is worth 
noting, however, that the hard-core condition (12. 3 p is 
also typically violated by density- functional theories |35| . 
The inability of approximation RFA to guarantee that 
gij (r) = for r < aij will be remedied by approximation 
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FIG. 1. (Color online) Plane A vs a^jox showing the shaded 
region — cr2/(o"i +0-2) < A < 2a"2/(o"i +0-2) where the first two 
singularities of gij{r), according to approximation RFA, are 
aij and Uik — a^j with k ^ j. The circles denote the systems 



>^^(^)-%-^^^^^.(^), 



(3.30) 



and Do{s) is the determinant of the matrix Q~^(s). Ex- 
plicit expressions of $ij(s) and Do{s) for binary mixtures 
arc given in Appendix [Cl 

Taking the Laplace inversion of Eq. (|3.27p . one finds 
that, in the interval < r < max(CTij, Tij) + e, 

H x^Q{r - nj)j,^j{r - Tij), (3.31) 

where 4>ij{r) and jikjir) are the inverse Laplace trans- 
forms of ^ij{s) and Tikj{s), respectively. 



Note that 



,(0) 



lim.. 



,^tj{s) 



L^]\ while 



analyzed in Sec. [V] 



7ifcj"(0) = lims->oo rifcj(s) = 0. Therefore, the contact 
values are 

2^(1) 2;r 

9iMt]) = -^ + a;„e(o-ij-Ty)7,«,j(o-y-Ty). (3.32) 

(jij (j^j 

As expected, Eq. p.32p reduces to Eq. (j2.40p in the ad- 
ditive case. 



RFA+ described in Sec. HITCl 



C. Approximation RFA^ 



3. Short-range behavior 

Before presenting approximation RFA-i-, we will need 
to restrict ourselves to cases where the first two singu- 
larities of gij{r), as given by approximation RFA, are Uij 
and Tij = min(crifc — akj\ k = 1, . . . ,n;k y^ j). As proven 
in Appendix |B1 the above requirement in the binary case 
{n = 2) implies the constraint —crilijJi -I- a-i) < A < 
2ail(a\ -f (T2), where, without loss of generality, it has 
been assumed 172 < <J\. This region of applicability is 
shown in Fig. [TJ 

Appendix [C] gives the expressions for gij (r) in the 
range < ?' < niax(crij, Ty) 4- e, where e is any positive 
value smaller than the separation between max((Ty- , Tij ) 
and the next singularity of gij(r), provided by approxi- 
mation RFA for binary mixtures. Extending to general 
n the arguments presented there, we can write 

Gy(s)=e-"-^$„-(s)H-2^pa;«e-"-*r,«j(s)-H--- , (3.27) 

where A: = k is the index corresponding to Tij, i.e., 
Tij = cTii^ — Qkj, and the ellipsis denotes terms headed 
by exponentials of the form e^"^* with A > ina.x{(Jij,Tij). 
In Eq. fiJTm . 



*y(*) = -^L^j{s)Qjj{s), 



r^fcj (s) 



1 L,k{s)Nkjis) 
s^ Do{s) ■ 



(3.28) 



(3.29) 



This new option for gij(r) will differ from approx- 
imation RFA only in the region min(CTij,Tij) < r < 
max(cry, Ty). More specifically, 

ffy(^)lRFA+ = 5»j('-)Irfa + — ^-'k Mr - <^tj) 

-8(r - Tjj)] 7i„j(r - Tij). (3.33) 

On account of Eq. p.3ip . Eq. p.33p can be equivalently 
rewritten as 



9t]ir)\ 



RFA4 



'^ir-<^^j) gz] l*^) I RFA ' '^y- < ^U" ' 

5y WIrfa + ©(»'- CTjj)0(t"jj - r) 

(3.34) 
We see from Eq. p.34p that the idea behind approxima- 
tion RFA+ is two-fold. On the one hand, it removes the 
unphysical violation of the property gij (r) = for r < aij 
that is present in option RFA when Tij < aij. On the 
other hand, if t^ > aij, approximation RFA-|_ extrap- 
olates to the region aij < r < Tij the functional form 
of gij (r) provided by approximation RFA in the region 
between Tij and the next singularity. 

In the interval < r < max(cry , Tij) + e. 



5yWI 



RFAi 



= -e(r-crij) [(/)jj (r-cr.y )+27rpxK7iKj {r 



In particular. 



(3.35) 



-(1) 



9v(<^ij)\RFA ^-^ + x^-f,^j{a,j-T,j). (3.36) 

+ (7ij (Tij 



D. Approximation RFA 



(m) 



T'ij > <Jij. In summary, option RFAl" is dcfinGd by 



^{r- cr^J ) g^J (r) \ j^p^ , r„- < cr„- , 



5»jWI 



RFAi 



5u(^)Irfa + ©(''- cr,j)e{nj - r) 

(3.37) 
Consequently, the contact values are 



In approximation RFA_|- the full functional form of 
likjir) is used. This can create some artificial problems 
in the region aij < r < nj when r^ > atj and the dis- 
tance Tij — fjij is rather large (as happens in the WR 
model). Reciprocally, if Ty — Uij is not large, it becomes 
unnecessarily complicated to consider the entire nonlin- 
ear function likjir) in the interval aij < r < Tij. Thus, 
we now propose a variant of approximation RFA+ , here 
denoted as RFA^"^ , whereby the full true function ^ji^j (r) 
is preserved if r^ < cr.y (in order to enforce the physical 
constraint of a vanishing RDF for r < aij) but is replaced 

by its mth degree polynomial approximation 7,-™ (?") if 



r(l) o 

Oy CTjj 



+ 



0(CTij -Ty)7iKj(cr.y 



,) 



-he(Ty - o,j)ii^>{a,j - r„)p.38) 

The polynomial y^^Hr) is obtained by truncating af- 
ter r"* the expansion of ^iujir) in powers of r. Such an 
expansion is directly related to that of the Laplace trans- 
form Vikjis) in powers of s~^ . For large s, Vi^jis) can 
be shown to be given by 



rzfcj (s) 
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(3.39) 



Consequently, the linear and quadratic approximations are 
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(3.40) 



tSW = 71,M+ i 



ifcj^ 



ifc 



r (0) ^fcj ^ (1) 
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bkj -L 
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L^bk, - lW 
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(3.41) 



Of course, the three sets of approximations RFA, (|3.36p . and p.38p in approximations RFA, RFA_|-, and 



RFA+, and RFA!;^'' reduce to the FY solution in the 
additive case. Obviously, RFA+ = RFA^°°'. In Sec. |V] 
we will generally use RFA\_ ' . 



IV. COMPARISON WITH MONTE CARLO 

SIMULATIONS FOR BINARY MIXTURES. THE 

EQUATION OF STATE 



RFA^ , respectively. 

In the compressibility route, the isothermal compress- 
ibility X is obtained from Eq. p. lip , where hij{0) = 



(1) 



p^/x^hij{0) = -Airp^/xlxJH^j , H-j' being the coef- 
ficient of s^ in the series expansion of s^Gij (s) in powers 
of s [cf. Eq. ((233)) ]. We recaU that Gij{s) is given by Eq. 
p.4p in approximation RFA. In approximations RFA-(_ 

and RFA^"\ Eqs. ((XMl) and ^^?ST\ imply that 



The compressibility factor Z is obtained via the virial 
and compressibility routes by Eqs. (|2.7p and (|2.13p . re- 
spectively. In the case of the virial route one needs the 
contact values of the RDF, which are given by Eqs. (|3.32p . 
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FIG. 2. (Color online) Compressibility factor as a function of 
the nonadditivity parameter for a symmetric binary mixture 
of NAHSs at pa"^ — 0.2 and two different compositions. The 
MC data are taken from Refs. [39, ISTIl . 
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(..) - ^^J 
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- 2TrpXi^ 



dri 



&{(7ij - Tij) 



(m). 



X7»Kj(r - r,j) + e(Ty - <7^J)r^^/{r - Ty) 



(4.2) 



In any case, for the sake of simplicity, we will restrict 
ourselves in most of this section to approximation RFA. 



A. Dependence of the EOS on nonadditivity 

Here we study the dependence of the EOS on the non- 
additivity parameter A by fixing all the other parameters 
of the mixture (density, composition, and size ratio). 



1. Symmetric binary mixtures 

Symmetric mixture are obtained when ai ~ <72 — cr. 
Therefore, in the additive case (A ~ 0) one recovers 
the one-component HS system, i.e., gii{r) ~ g22(r) ~ 
gi2i'r) = g{r), regardless of the value of xi. 

Figure [5] compares the compressibility factor obtained 
from MC simulations [3a, |33| with that predicted by 
approximation RFA for some representative symmetric 
systems. We observe that approximation RFA repro- 
duces quite well the exact simulation data at all val- 
ues of the nonadditivity parameter. At this low den- 
sity {pa^ = 0.2,77 ~ 0.105) the virial and compressibility 
routes are not distinguishable on the scale of the graph. 



2. Asymmetric binary mixtures 

Asymmetric mixtures correspond to cti 7^ (T2- In that 
case, when A = one recovers the AHS mixture. 



FIG. 3. (Color online) Compressibility factor as a function 
of the nonadditivity parameter for an equimolar asymmetric 
binary mixture of NAHSs with a size ratio a^jcn = 1/3 at a 
packing fraction r\ = 0.5. The symbols [v] and [c] stand for 
the virial and compressibility routes, respectively. The MC 
data are taken from Ref. [38l |. 



Figure [3] shows the A dependence of Z for negative 
nonadditivity and an equimolar {x\ ~ X2 = 5) ^-sym- 
metric mixture {cr 2/(71 = 1/3) at a relatively large den- 
sity (77 = 0.5). In this case the virial route of approx- 
imation RFA underestimates the values of Z ^ while the 
compressibility route overestimates them. This is also a 
typical behavior of the FY equation for AHSs. It is thus 



tempti ng to try the Z = ^Z"" -\-^Z'^ interpolation recipe 



[39|-|42[ , which is known to work well in the additive case. 
From Fig.|3]we see that indeed the interpolation formula, 
as applied to approximation RFA, reproduces quite well 
the exact simulation data, except for A < —0.8. 



B. Dependence of the EOS on the size ratio 

Next, we study the dependence of Z on the size ratio 
(J2I o\ by fixing all the other parameters of the mixture 
(density, composition, and nonadditivity). 

The three panels of Fig. 0] show Z vs (T2/o'i for a 
slightly negative nonadditivity A = —0.05 (top panel), a 
moderate positive nonadditivity A = 0.2 (middle panel), 
and a larger positive nonadditivity A = 0.5 (bottom 
panel). We observe again that the interpolation recipe 
Z = ^Z"" + ^Z'^ for approximation RFA agrees well with 
the exact simulation data, with the exception of a region 
close to the size symmetric mixture {02! <J\ — V) for posi- 
tive nonadditivity and moderate density (middle panel). 



C. Contact values 

In Sec.|V]we will analyze the RDF gijir) predicted by 
approximations RFA and RFA!^ . Before doing so, and 
as a bridge between the thermodynamic and structural 
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FIG. 4. (Color online) Compressibility factor as a function of 
the size ratio a2/(Ji for binary asymmetric NAHS mixtures 
with a;2 = |, A = —0.05, and rj = 0.5 (top panel); X2 = j, ^, 
A = 0.2, and rj = 0.2 (middle panel); X2 = j,^, A = 0.5, 
and ?7 = 0.075 (bottom panel). In the bottom panel only 
the theoretical data obtained from the virial route are shown 
since they practically coincide with those obtained from the 
compressibility route. The MC data are taken from Ref. [3g|. 



TABLE I. Contact values for some binary equimolar symmet- 
ric NAHS mixtures. The MC and PY data were taken from 
Ref. [J]. The labels correspond to systems common to those 
listed in Table HH 



Label 



pa Source giiji^ ) gi2{o'i 



D 



0.05 0.8 



0.0 



0.8 



-0.05 0.8 



A 



B 



-0.1 



-0.3 



-0.5 



1.0 



1.0 



1.0 



MC 


5.305 


3.762 


PY 


4.451 


3.516 


RFA 


4.006 


3.617 


rfaV'' 


4.580 


3.617 


MC 


3.971 


3.971 


PY 


3.581 


3.581 


RFA 


3.581 


3.581 


rfaV^' 


3.581 


3.581 


MC 


3.117 


3.801 


PY 


2.925 


3.394 


RFA 


2.971 


3.148 


rfa(!' 


2.971 


3.445 


MC 


3.394 


5.363 


PY 


3.209 


4.395 


RFA 


3.497 


3.883 


rfaV^' 


3.497 


4.763 


MC 


2.168 


2.798 


PY 


2.141 


2.543 


RFA 


2.441 


2.251 


rfa(!' 


2.441 


2.875 


MC 


2.103 


1.528 


PY 


2.060 


1.493 


RFA 


2.139 


1.407 


rfaV^' 


2.139 


1.279 



biliary symmetric mixtures th = T22 = 1712 = o'(l + A) 
and T12 = (7, it turns out that gii(cr+) = 522(0"^) is com- 
mon in approximations RFA and RFAI^'^ if A < 0, while 
512(0" 12) is common in both approximations if A > 0. 

From Table |T] wc observe that approximation RFA^ ' 
is superior to the PY theory in estimating the true con- 
tact values, both for positive and negative nonadditivity, 
except in the cases of gii{a~^) for pa^ = 1 and A = —0.3 
and of gi2{(7t2) foi' P'^^ = 1 s-i^d ^ = —0.5. 



V. COMPARISON WITH MONTE CARLO 

SIMULATIONS FOR BINARY MIXTURES. THE 

STRUCTURE 



properties, it is worth considering the contact values. Ta- 
ble U provides the contact values for some binary equimo- 
lar symmetric NAHS mixtures (cti = a2 = cr, xi = X2 ^ 
i), as obtained from MC simulations [J, numerical solu- 
tions of the PY integral equation [j| , and our approxima- 
tions RFA [Eq. (13311)] and RFA^^^ [Eq. (^351) ]. Since for 



The RDF of approximation RFA is analytically and 
expHcitly given in Laplace space by Eqs. p.4p - p.6p and 
p.l5p - p.2ip . In real space, rg,y(r) is easily found by 
taking the inverse Laplace transform of Gij{s) through 
the numerical scheme described in Ref. [4J]. To get gij{r) 

in approximation RFAl^ , one needs to make use of Eq. 

(|3.37p . where 7|fcJ (j^) is explicitly given by Eqs. p.40|) 
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TABLE II. The six binary NAHS mixtures considered in the 
analysis of the structure. The last column gives the compress- 
ibility factor as obtained from our MC simulations. 



Label o^jox 



xi 



pal 



A 


1 


-0.1 


1/2 


1.0 


0.5236 


8.648 


B 


1 


-0.5 


1/2 


1.0 


0.5236 


3.429 


C 


4/5 


-0.444 


1/3 


1.0 


0.3533 


2.335 


D 


1 


0.05 


1/2 


0.8 


0.4189 


9.083 


E 


1 


0.25 


1/2 


0.3 


0.1571 


2.556 


F 


4/5 


0.25 


1/3 


0.3 


0.1060 


1.876 




FIG. 5. (Color online) RDF for system A of table |TI1 



and p.4ip for m — 1 and to = 2, respectively [45|. No- 
tice that, while the true RDF has to be symmetric under 
exchange of species indices, the RDF obtained from ap- 
proximation RFA or RFA+ is, except for symmetric and 
equimolar mixtures, not symmetric, i.e., gij{r) ^ gjiir) 
if i 7^ J. Although this artificial asymmetry is gener- 
ally small from a practical point of view, it represents 
a penalty we pay for our extension of the AHS solution 
of the PY equation. To cope with this shortcoming, we 
just redefine the like-unlikc RDF as the symmetrized one 
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In a binary mixture. Til — fi2+ai2 — o'i + i(o'i+cr2)A, 125/128 ~ 0.98. 
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FIG. 6. (Color online) RDF for system B of table [ 



T22 = cri2 - ai2 = (T2 + ^{o-i + cr2)A, and ri2 = ^(cti -I- 
a2)- Therefore, th < cti and r22 < cr2 for A < 0, while 
^12 < cri2 for A > 0. In what follows, we will truncate 

In order to evaluate the merits and limitations of the 
structural properties predicted by our approximations, 
we have performed canonical MC simulations of the bi- 
nary NAHS system with N = 2196 particles and lO^N 
MC steps per run. The cell index method has been used 
\4&j . The statistical error on the RDF is within the size 
of the symbols used in the graphs reported. 

We have chosen six representative systems, all within 
the region — a2/(o'i + a2) < A < 2tT2/(ci -I- CT2) assumed 
in the construction of approximation RFA_|.. Those six 
systems are represented in Fig.[l]and their respective val- 
ues of composition and density arc displayed in Table |lll 
Three of the mixtures have a negative nonadditivity (A, 
B, and C), while the other three have a positive nonaddi- 
tivity (D, E, and F). Moreover, there are four equimolar 
symmetric mixtures (A, B, D, and E) and two asym- 
metric ones (C and F). In those two latter cases, how- 
ever, both species contribute almost equally to the (nom- 
inal) packing fraction 77 since a;icri/a:2cr2 ~ (5/4)'^/2 = 
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FIG. 8. (Color online) RDF for system D of table HI 



FIG. 7. (Color online) RDF for system C of table [U 



gi2(r) in the interval cti2 = 0.9cti < r < T12 = cri. In fact, 
approximation RFA presents an artificial discontinuity of 
the first derivative g'12 ('') at r = cri . This is corrected by 
approximation RFAl|_ , which presents a good agreement 
with the MC results for r < ai. In spite of this, we 
observe that approximation RFAI^ underestimates the 
contact value 312 (('■12)1 in agreement with the entry of 
Table H] corresponding to case A. 



A. Negative nonadditivity 



1. Symmetric m,ixtures 

Figures [5] and [6] display the RDF for systems A and B, 
respectively. System A is only slightly nonadditive and 
we observe that both approximations RFA and RFA^ 
do a very good job. On the other hand, while RFA and 
RFA^' coincide for gii{r) with r > ai, they differ for 



In the case of system B the nonadditivity is larger and, 
according to Fig. [HI the performance of our approxima- 
tions is still good for .gii(r) but worsens for gi2{r). In 
fact, gi2('')|RFA turns out to be better than 9i2(j')\„pj^(i) 
in the region ai2 ~ O.Scri < r < T12 = cti, in agreement 
with the entry of Table U corresponding to case B. In 
any case, it is interesting to remark that approximation 
RFA^ succeeds in capturing the non-monotonic behav- 
ior of 5i2(?') very near r ~ ai2 observed in the simula- 
tions. 
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FIG. 9. (Color online) RDF for system E of table [ 



2. Asymmetric mixture 

The only case representing an asymmetric mixture 
with negative nonadditivity (system C) is shown in Fig. 
[71 Again, the MC like-like RDF are very well reproduced 
by the two approximations. In the case of the like-unlike 
function gi2{r), approximation RFA^ clearly improves 
approximation RFA in the region cri2 ~ O.Scri < r < 
Ti2 = 0.9cri. Apart from that, both approximations over- 
estimate gi2{i') between r = T12 = 0.9cri and the location 
of the first minimum at about r ~ 1.25o'i. In Fig. [7] we 
have taken gi2{r) — >■ \[gi2{r) + g2i{i')], as explained at 
the beginning of this section. Prior to this symmetriza- 
tion, the maximum relative deviation between gi2 (r) and 
g2\{f') occurs at r ~ 0.75(Ti and is less than 5%. 



B. Positive nonadditivity 

1. Symmetric mixtures 

Let us consider now positive nonadditivities, start- 
ing with symmetric mixtures. Figures [5] and [5] show 
the results for systems D and E, respectively. For a 
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FIG. 10. (Color online) RDF for system F of table HH 



small nonadditivity A = 0.05, both approximations pro- 
vide very good results, except for gii{r) near contact 
(see also Table |l|. Notice, however, that approximation 

RFA^ improves approximation RFA in the narrow re- 
gion (Tl < 7' < Til = 1.05(Tl. 

For a larger nonadditivity (system E), Fig.|9]shows the 
excellent job made by approximation RFA^_ in the inter- 
val CTi < r < Til = 1.25(71. In the case of the like-unlike 
correlation function, however, the approximations over- 
estimate the values between tTi2 and the first minimum 
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FIG. 11. (Color online) RDF for the WR model at pafi 
0.28748. The MC data are taken from Ref. H. 



FIG. 12. (Color online) RDF for the WR model at pcrfa = 0.4. 
The MC data are taken from Ref. El. 



(r ~ 2ai). 



The Widom— Rowlinson model 



2. Asymmetric mixture 

Figure [10] displays the three functions gij (r) for the 
asymmetric system F. As in case E, approximation 

RFA^ nicely reproduces the exact results from the sim- 
ulation for the like-like correlations and corrects the un- 
physical kink of approximation RFA occurring at rn ~ 
1.2250-1 and T22 = 1.025tTi. Interestingly enough, al- 
though the values of A and paf are the same in systems E 
and F. the performance of the approximations for (712 (j") 
is much better in case F (asymmetric mixture) than in 
case E (symmetric mixture). This might be partially due 
to the fact that the packing fraction 77 is smaller in sys- 
tem F than in system E. For the asymmetric system F, 
we have found that the maximum relative deviation be- 



tween gi2(r) and 52i('') takes place at r = (T12 
is less than 0.5%. 



jCTi and 



As recalled in Sec. HI the WR model corresponds to 
an equimolar symmetric binary NAHS mixture where 
ci = 172 = and (T12 j^ 0. The model is then fully 
characterized by the reduced density, pc7i2- The critical 
demixing reduced density for this model is around 0.75 

MM- 



The nonadditivity parameter of the WR model is A = 
add 1 — >. 00, so it lies well outside the "safe" region 



^12/^^2 



for our approximation RFA+ (see Fig. [I}. To compen- 
sate for this, we replace here approximation RFA\ by 

(2) 
approximation RFA^ . 

We see from Figs. [TT] and [1^ that approximation 

(2) 

RFA^^ does a much better job than expected at the 
two densities considered. The main drawbacks of the 
theory are that the contact value gii(O) is dramatically 
overestimated and the behavior of gi2{r) for r > CT12 is 

qualitatively wrong. In spite of this, it is remarkable that 

(2) 
approximation RFAY captures well the global properties 

of the RDF in this extreme system. 
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VI. SUMMARY AND CONCLUSIONS 

The importance of the NAHS model in hquid state the- 
ory cannot be overemphasized. When the reference or 
effective interaction among the microscopic components 
(at an atomic or a colloidal level of description) of a sta- 
tistical system is modeled as of hard-core type, there is 
no reason to expect that the interaction range aij corre- 
sponding to the pair (i, j) is enslaved to be the arithmetic 
mean of the interaction ranges ct^ and CTj corresponding 
to the pairs (i, i) and (j, j), respectively. Therefore, in an 
n-component NAHS mixture the number of independent 
interaction ranges is n{n -\- 1)/2, in contrast to the num- 
ber n in an AHS mixture. It is then not surprising that, 
while an exact solution of the PY theory exists for AHS 
systems |13| , numerical methods arc needed when solving 
the PY and other integral-equation theories for NAHSs 
[J|. Therefore, analytical approaches to the problem can 
represent attractive and welcome contributions. 

In this paper we have constructed a non-perturbative 
fully analytical approximation for the Laplace transforms 
Gij{s) of rgij(r), where gij{r) is the set of RDF of a gen- 
eral 3D NAHS fluid mixture. Our approach follows sev- 
eral stages. The starting point is the analytical PY solu- 
tion for AHSs, Eqs. ([236)) - ([238)) . Exploiting the connec- 
tion between the exact solutions for ID NAHS and AHS 
mixtures [see Eqs. ([^T5)) and ^^], the AHS PY solu- 
tion is rewritten in an alternative form, Eqs. p.ip - p.4p . 
Our approximation RFA consists of keeping the form 
p.4p . except that af^^ in Eq. p.ip is replaced by ai, [cf. 
Eq. ()3.5|) ] and ai in Eq. p.3p is replaced by bij = Uij 



'I] 

cm +ai 



(0) 



Moreover, the parameters Li-' and L, 



(1) 



[cf. Eq. 

are no longer given by Eq. ()2.4ip but arc determined by 
enforcing the condition ()2.4p or, cquivalently, Eq. ()2.33p . 
This results in Eqs. p.l5p - p.2ip . and so the problem be- 
comes completely closed and analytical in Laplace space. 
The equation of state is obtained either via the virial 
route ()2.7p through the contact values ()3.32p or via the 
compressibility route ()2.1ip through the coefficients H^- 
in the expansion of s'^Gij(s) in powers of s, Eq. ()2.33p . 

The penalty we pay for "stretching" the AHS PY so- 
lution to the NAHS domain in the way described above 
is that gij (r) may not be strictly zero for r < cr,;j or may 
exhibit first-order discontinuities at artificial distances. 
To deal with this problem, we have restricted ourselves 
to mixtures such that the first two singularities of gij (r) 
are (Tij and t^ = n\m{f7ik — akj ;k — 1, . . . ,n;k ^ j). In 
the binary case {n ~ 2) this restriction corresponds to 
-cr2/(cri -h (72) < A < 2cr2/(cri + 02) (scc Fig.[T|). Next, 
we have constructed a modified approximation RFA+ 
whereby either gij (r) is truncated for r < aij if Tij < Uj 
or the behavior of gij (r) for r > Tij is extrapolated to the 
interval Gij < r < Tij if Tij > aj [cf. Eq. ()3.34p ]. From 
a practical point of view, the latter extrapolation can 
be replaced by a polynomial approximation (e.g., linear 
or quadratic), yielding approximation RFA_|^" [cf. Eq. 
p.37p ]. This is sufficient to guarantee that the slope of 



gij{r) is continuous everywhere for r > iTij. 

For comparison with MC data of the equation of state 
we have used approximation RFA since its local limita- 
tions at the level of the RDF are largely smoothed out 
when focusing on the thermodynamic properties. The 
results show that, if the density is low enough as to make 
both thermodynamic routes practically coincide, our ap- 
proximation accurately predicts the MC data, as shown 
in Fig. [2] and in the bottom panel of Fig. |4| For larger 
densities, the virial and compressibility routes tend to 
underestimate and overestimate, respectively, the simu- 
lation values, this being a typical PY feature. As in the 
AHS case, the simple interpolation rule Z = ^Z"" + ^Z'^ 
provides very good results, except for large nonadditiv- 
itics (see Fig. |3| and the top and middle panels of Fig. 

Regarding the structural properties, approximation 
RFA^_ ' is found to perform quite well. The contact values 
are generally more accurate than those obtained from the 
numerical solution of the PY integral equation, at least 
for symmetric mixtures, as shown in Table |I| Compar- 
ison with our own MC simulations shows a very good 
agreement, except in the case of the like- unlike RDF for 
distances smaller than the location of the first minimum 
for large nonadditivities (cf. Figs. ISI ITU)) . On the other 
hand, even in the case of the WR model (A — > 00, well 

beyond the "safe" region of Fig. |T|) our approximation 

(2) 
RFAY does a much better job than expected, as illus- 
trated in Figs. [TT] and [T^ 

In conclusion, one can reasonably argue that our ap- 
proximation RFA, along with its variants RFA-|_ and 
RFA^" , represent excellent compromises between sim- 
plicity and accuracy. We have tried other alternative 
analytical approaches (simpler as well as more complex) 
also based on the PY solution for AHSs but none of them 
has been found to present a behavior as sound and consis- 
tent as those proposed in this paper. We expect that they 
can be useful in the investigation of such an important 
statistical-mechanical system (both by itself and also as 
a reference to other systems) as the NAHS mixture. 

The work presented in this paper can be continued 
along several lines. In particular, we plan to explore in 
the near future the predictions for the demixing transi- 
tion from our approximations. It is also worth exploring 
the NAHS theory that arises when the starting point is 
not the PY solution for AHSs but the more advanced 
RFA proposed in Ref. [ij], which contains free parame- 
ters that can be accommodated to fit any desired EOS 
in a thermodynamically consistent way. 
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Appendix A: Explicit expressions of Gij{s) for binary 
mixtures in approximation RFA 



By performing the inversion of the matrix p.8p and 
carrying out the matrix product in Eq. p.4p one gets 



Gii(s) 



Dis) 



Luis) 



2TTpX2 , ■ 

1 :T—N22[s) 



+ ^Lu(.)i22(.s)e-(-+-)^ 



Li2(s)L2i(s)e-2'^^^-^+^^Li2(s)A^2i(s)e-(^-+'^-)-' 



(Al) 



Gl2(s) 



D{s) 



Li2{s) 



1 - 



2Trpx 



-Nn{s) 



e---^ + ?^Ln(s)iVi2(s)e-(-^+-^)^/2 \ , 



(A2) 



where the quadratic functions Nkj{s) can be found in Eq. p.9p and 

{2T:pyxiX2 



Dis) = 



1 '^^P^^ AT ( \ 

1 ^— Afii(.s) 



27rpx2 

1 ^— A^22(s) 



Nl2{s)N2l{s) 



2-Kpxx . , 
-— ^iii(s) 



1- 



2'KpXi 



N22is) 



.--+"=^L22is) 



2'KpXl 



Niiis) 



4:TT pXlX2 



Lii{s)L22{s)e 



-~{<yi+a-2)s 



Li2{s)L2i{s)e 



— 2ai2S 



+ii2(s)A^2i(s)e-("^^+''^^)^ + i2i(s)A^i2(.s)e-("^^-''^^)^ 



(A3) 



is the determinant of the matrix Q ^. The expressions 
for 6*22(5) and G2i{s) can be obtained by the exchange 
1 o2. 



Appendix B: Ordering of singular distances in 
approximation RFA for binary mixtures 

By "singular" distances we will refer to those values 
of r where the RDF gij (r) or any of its derivatives have 
a discontinuity. Physical singularities are located, for 
instance, at r = cr.y and r = Oik + (Xfej, k = l,...,n. 
Apart from that, approximation RFA introduces spurious 
singularities at other distances. 

Let us particularize to a binary mixture. The physical 
leading singularity of gij (r) should be located at r = cTij . 
However, according to Eq. (|A1[) . the leading singularity of 
gii{r) takes place at r = min(cri, ai2 + 012, 2ai2). Analo- 
gously, the leading singularity of 1722 {t) is located at r = 
min(cr2,(Ti2 — ai2, 2cri2). Finally, Eq. (jA2[) shows that the 
leading singularity of gi2{'r) is r = ■imin(2cri2,(Ji + (72). 
Note that we have assumed cri2 — ai2 > 0, so that the 
denominator D(s), Eq. (|A3p . does not affect the leading 
singularity of gij{r). 

It is thus important to determine the relative order- 
ing of the values cti, 02, cti2 — ai2! cri2 -I- ai2, 2ai2, and 
ai+(T2. Such an ordering depends on the values of A and 



I 

R = a'2/ci, where, without loss of generality, we assume 
that (72 < CTi. A detailed analysis shows that the A-R 
plane can be split into 13 disjoint regions with distinct or- 
der for the above singular distances. Those regions arc in- 
dicated in Fig.[T21 while Table Hill shows the order apply- 
ing within each region. Note that ai2 — a\2 is negative in 
Regions Ilf, Ilg, and Ilh, i.e., if -1 < A < -2i?/(l -^ i?), 
thus invalidating those regions from the preceding anal- 
ysis. 

We observe that a\ and 02 are indeed the leading sin- 
gularities of (711 (r) and 322 (j"), respectively, for positive 
nonadditivity (regions la-Ie). Reciprocally, (j\2 is the 
leading singularity of 5i2('') for negative nonadditivity 
(regions Ila-IIh). 

In order to construct approximation RFA-|_, we want 
to restrict ourselves to those regions such that the two 
leading singularities of gwir) are a\ and rn = a\2 + 012. 
Inspection of Table IIIII shows that Regions Ilc-IIh are 
discarded by this criterion. In the remaining regions 
the leading singularity of .911(7") is min((Ti,cri2 -f 012) 
but the next one is not necessarily max((7i,cri2 + 012) 
since the latter value competes with min((Ti, 1712 + 012) + 
min(cr2,(Ti2 — ai2,2cri2), where the term min(o'2,o'i2 — 
012,2(712) comes from the denominator D[^s) [cf. Eq. 
(|A3p ]. It can be checked that max((7i,(7i2 + 012) > 
min((7i,t7i2 + 012) -I- min(cr2,o'i2 — ai2,2cri2) in Regions 
Ic-Ie. Therefore the two first singularities of gwir) are 
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TABLE III. Order of the singular distances cri, (T2, cri2 — ai2, (T12 + ai2, 2(7i2, and ai + (J2 in each of the regions of Fig. 1131 



Region 



Order 



la 

lb 

Ic 

Id 

le 

Ila 

lib 

lie 

lid 

He 

Ilf 

Ilg 

Ilh 



0"12 
0"12 
0-12 














- ai2 

- (112 

- CLV2 



< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



02 
02 
02 
02 
02 



''"12 
(7l2 
0"12 



- 0,12 

- ai2 

- ai2 



o"i2 — ai2 
(T12 — ai2 









< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



(^12 — "12 

o"i 

0"12 — ai2 

o\ 

02 

0'12 + ai2 

C2 
C"12 + ai2 

(T12 + ai2 
02 

2cri2 

2cri2 



< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



a\2 — a\2 
(y\ 

0"12 — ai2 

(Jl + 0-2 

(y\2 + "12 

0"2 
0"12 + ai2 

02 

2cri2 

2cri2 

0-2 

''"12 + "12 



< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



0"12 + ai2 
0"12 + ai2 

fTl + (72 
(Tl + (T2 

0"12 — ai2 
o"i 

O"! 

2cri2 

2cri2 

(72 

(712 + ai2 

(712 + ai2 

(72 



< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



(7l + (72 
(71 + (72 

(712 + ai2 
(712 + ai2 
(712 + ai2 

2(712 
2(712 

(71 

(71 

(71 

(71 

(71 

(7l 



< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 
< 



2(712 

2(712 

2(712 

2(712 

2(712 

(7l + (72 

(7l + (72 

(7l + (72 

(71 + (72 

(71 + (72 

(71 + (72 

(71 + (72 

(7l + (72 




■1.0 



Ila -1{\-Ry{\+K) 
-RI{\^K) X lib 



.0 5 \ -27?/(l+i?) 
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-1/(1+7?) 
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FIG. 13. (Color online) Plane A vs _R = (72/(71 showing the 
regions with different ordering of the distances (7i , (72 , (J\2 — 
ai2, (712 + ai2, 2(712, and (7i + (72. 



tTi and Til = cri2 + 0,12 in Regions la, lb, Ila, and lib 
only. It turns out that in those four regions the two lead- 
ing singularities of g22(j') are CT2 and T22 = cri2 — ai2, 
and the two leading singularities of §12(1") are (T12 and 

T12 = 2(^1 +0'2)- 

In summary, Regions la, lb, Ila, and lib are the only 
ones where the two leading singularities of gij (r) are ct^ 
and Tij = (Tjfc - flfcj with k ^ j. 



Appendix C: Short-range forms of gij (r) for binary 
mixtures in approximation RFA 

In what follows we assume that — (T2/(cri + (72) < A < 
2cr2/(cri + (72), which corresponds to Regions la, lb, Ila, 
and lib of Fig. [T31 As discussed in Appendix [B] this 
guarantees that the first two singularities oi gij{r) arc (Ty 
and Tij = Uik — cikj with k ^ j . The aim of this Appendix 
is to give the expressions of gij{r) in the region < r < 
max ((Tij, Tij) + e, where e is smaller than the separation 
between max(CTy, Ty) and the next singularity. 

It is convenient to assign a bookkeeping parameter z to 
e^"; so that, for instance, e~"''^'^ becomes z°''^e~"''^'' . We 
will set 2 = 1 at the end of the calculations. Therefore, 
the denominator D{s) given by Eq. (jA3p becomes 



D{s) = Do{s)+o(z"), 



(CI) 



where 
Do{s) 



{27:p)^xiX2 



Nn{s) 



2'KpX2 



N22{s) 



Nl2{s)N2l{s). 



(C2) 



In Eq. (jCip . o(z'"') denotes terms that are negligi- 
ble versus z" in the (formal) limit z — > 0, i.e., 
lim^-^o z~"o(z") = 0. From Eq. (|A1[) we see that the two 
leading terms in Gii(s) are of orders z'^^ and z'^^^+'^i^: 

Gii(s) = $n(s)e-"i^z'^i + 27rpx2Ti2i{s)e-^^'''^^''^' 



xz 



""12+112 



+ o(z'"i) + o(z 



''"12+112 



), 



where 

$ii(s) 



Dois) 



Luis) 



1 ^— /V22(s) 



(C3) 



(C4) 



and rifej(s) is given by Eq. p.29p . Analogously, 
Gi2(s) = $i2(s)e-"i=^z"i^ +27rpxirn2(s)e-("^+"^)^/' 

X2('7l+<72)/2 _,_ ^(^<T12) _^ o(2:('^l+'-2)/2)^ ((.5) 
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G2i(.s) = <p2iis)e-''"'z'''' + 27rpx2r22i(s)e-('^i+"^)^/2 



G22(s) = $22(s)e-"^^z"^ +27rpxir2i2(s)e-("i^-''")^ 
where 



*12(s) 



^o(s) 



il2(s) 



2'Kpx 



-Nn{s) 



(C8) 



522(?') = -6(r - cr2)022(r - (72) H 6(r - (T12 + ai2) 

r r 

X7212(^-cri2 + 012), (C14) 

where we have already set z = 1. In Eqs. (jClip - ()C14p . 
4>ij{r) and "fikjir) are the inverse Laplace transforms of 
^ij{s) and rifej(s), respectively. 

Since 4'ij{0) — linis_i.oo ^tjis) = ^ij ^ ^^^ contact val- 
ues in approximation RFA are 



r(l) 



(C15) 



*2l(s) 



$22(5) 



Doi-s 



■L2iis) 



Dois) 



^22(5) 



2TTpX2 , ■ 

1 5— A^22(s) 



27rpxi ' 

1 ^— iVii(s) 



(C9) 



(CIO) 



Laplace inversion of Eqs. (|C3p and (|C5[) - (|C7p shows 
that in the interval < r < max{aij , r^ ) + e we obtain 



1 ^TT OX'^ 

gii{r) = -&ir - (Ti)(pii{r - cti) H — e(r - 0-12 - 012) 

r r 



'Xli2iir - <Ji2 -012), 



(Cll) 



3i2(?') = -B(r - (Ti2)(?>i2(r - 0-12) H e(r 



0-1 + cr2 , 
X7ii2(?- ^ ), 



^1 + ^2 X 

2 '^ 
(C12) 



512 K2) = ^^ + ^^ 

1712 0'12 



52i(.c^i2) ! 

0'12 0'12 



7112 I 0-12 ^ ^ ) , (C16) 



7221 I 0-12 ^ „ 1 , (C17) 



522(0-^) 



L, 



(1) 

22 



0-2 



(CIS) 



in Regions la and lb (A > 0). On the other hand, in 
Regions Ila and lib (A < 0), 



9iii<7t) = ^ 



l[^ 2ttpx2 



CTi CTi 



7121 (cti - Cri2 - 012) , (C19) 



ffl2(o-]*2) 



r(l) 
-^12 

CT12 



52l(a+2) = ^, 
0'12 



(C20) 



(C21) 



52i('') = -^{r ~ (T2i)4>2i{r - (T12) H B(r- ) 

r r ^ 



O"! + cr2 , 
X722l(?' -Z ), 



2 
(C13) 



L. 



922{a^) = ^^ + 



2TTpXi 



7212 (ct2 - cri2 + 012) . (C22) 



0'2 ^2 

A more compact form is provided by Eq. p.32p . 
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